Original Authors: Belinda Phipson, Anna Trigos, Matt Ritchie, Maria Doyle, Harriet Dashnow, Charity Law Based on the course RNAseq analysis in R delivered on May 11/12th 2016
Before starting this section, we will make sure we have all the relevant objects from the Differential Expression analysis present.
suppressPackageStartupMessages(library(edgeR))
load("Robjects/DE.Rdata")
We have a list of significantly differentially expressed genes, but the only annotation we can see is the Entrez Gene ID, which is not very informative.
edgeR provides a function plotSmear that allows us to visualise the results of a DE analysis. In a similar manner to the MA-plot for microarray data, this plot shows the log-fold change against log-counts per million, with DE genes highlighted:
summary(de <- decideTestsDGE(lrt.BvsL))
[,1]
-1 2957
0 11232
1 1615
detags <- rownames(dgeObj)[as.logical(de)]
plotSmear(lrt.BvsL, de.tags=detags)
However, on such a plot it would be nice to add labels to highlight the genes with most evidence for being DE, or our favourite genes. To perform such a task we need to map between the identifiers we have in the edgeR output and more familiar names.
Finally, we will look at sophisticated visualisations that allow us to incorporate information about the structure of a gene, level of sequencing coverage.
There are a number of ways to add annotation, but we will demonstrate how to do this using the org.Mm.eg.db package. This package is one of several organism-level packages which are re-built every 6 months. These packages are listed on the annotation section of the Bioconductor, and are installed in the same way as regular Bioconductor packages. An alternative approach is to use biomaRt, an interface to the BioMart resource. BioMart is much more comprehensive, but the organism package fit better into the Bioconductor workflow.
source("http://www.bioconductor.org/biocLite.R")
biocLite("org.Mm.eg.db")
# For Human
biocLite("org.Hs.eg.db")
The packages are larger in size that Bioconductor software pacakges, but essentially they are databases that can be used to make offline queries.
library(org.Mm.eg.db)
First we need to decide what information we want. In order to see what we can extract we can run the columns function on the annotation database.
columns(org.Mm.eg.db)
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID" "ENZYME"
[8] "EVIDENCE" "EVIDENCEALL" "GENENAME" "GO" "GOALL" "IPI" "MGI"
[15] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID" "PROSITE" "REFSEQ"
[22] "SYMBOL" "UNIGENE" "UNIPROT"
We are going to filter the database by a key or set of keys in order to extract the information we want. Valid names for the key can be retrieved with the keytypes function.
keytypes(org.Mm.eg.db)
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID" "ENZYME"
[8] "EVIDENCE" "EVIDENCEALL" "GENENAME" "GO" "GOALL" "IPI" "MGI"
[15] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID" "PROSITE" "REFSEQ"
[22] "SYMBOL" "UNIGENE" "UNIPROT"
We should see ENTREZID, which is the type of key we are going to use in this case. If we are unsure what values are acceptable for the key, we can check what keys are valid with keys
keys(org.Mm.eg.db, keytype="ENTREZID")[1:10]
[1] "11287" "11298" "11302" "11303" "11304" "11305" "11306" "11307" "11308" "11350"
It is a useful sanity check to make sure that the keys you want to use are all valid. We could use %in% in this case.
## Build up the query step-by-step
my.keys <- c("50916", "110308","12293")
my.keys %in% keys(org.Mm.eg.db, keytype="ENTREZID")
[1] TRUE TRUE TRUE
all(my.keys %in% keys(org.Mm.eg.db, keytype="ENTREZID"))
[1] TRUE
Let’s build up the query step by step.
## to be filled-in interactively during the class.
select(org.Mm.eg.db,
To annotate our results, we definitely want gene symbols and perhaps the full gene name. Let’s build up our annotation information in a separate data frame using the select function.
ann <- select(org.Mm.eg.db,keys=rownames(results),columns=c("ENTREZID","SYMBOL","GENENAME"))
'select()' returned 1:1 mapping between keys and columns
# Have a look at the annotation
ann
Let’s double check that the ENTREZID column matches exactly to our results rownames.
table(ann$ENTREZID==rownames(results))
TRUE
15804
We can bind in the annotation information to the results data frame. (Please note that if the select function returns a 1:many mapping then you can’t just append the annotation to the fit object.)
results.annotated <- cbind(results, ann)
results.annotated
We can save the results table using the write.csv function, which writes the results out to a csv file that you can open in excel.
write.csv(results.annotated,file="B.PregVsLacResults.csv",row.names=FALSE)
A note about deciding how many genes are significant: In order to decide which genes are differentially expressed, we usually take a cut-off of 0.05 on the adjusted p-value, NOT the raw p-value. This is because we are testing more than 15000 genes, and the chances of finding differentially expressed genes is very high when you do that many tests. Hence we need to control the false discovery rate, which is the adjusted p-value column in the results table. What this means is that if 100 genes are significant at a 5% false discovery rate, we are willing to accept that 5 will be false positives. Note that the decideTests function displays significant genes at 5% FDR.
Challenge
Re-visit the
plotSmearplot from above and use thetextfunction to add labels for the names of the top 200 most DE genes
Another common visualisation is the volcano plot which display a measure of significance on the y-axis and fold-change on the x-axis.
signif <- -log10(results.annotated$FDR)
plot(results.annotated$logFC,signif,pch=16)
points(results.annotated[detags,"logFC"],-log10(results.annotated[detags,"FDR"]),pch=16,col="red")
Before following up on the DE genes with further lab work, a recommended sanity check is to have a look at the expression levels of the individual samples for the genes of interest. We can quickly look at grouped expression using stripchart. We can use the normalised log expression values in the dgeCounts object (dgeCounts$counts).
library(RColorBrewer)
par(mfrow=c(1,3))
normCounts <- dgeObj$counts
# Let's look at the first gene in the topTable, Krt5, which has a rowname 50916
stripchart(normCounts["110308",]~group)
# This plot is ugly, let's make it better
stripchart(normCounts["110308",]~group,vertical=TRUE,las=2,cex.axis=0.8,pch=16,col=1:6,method="jitter")
# Let's use nicer colours
nice.col <- brewer.pal(6,name="Dark2")
stripchart(normCounts["110308",]~group,vertical=TRUE,las=2,cex.axis=0.8,pch=16,cex=1.3,col=nice.col,method="jitter",ylab="Normalised log2 expression",main=" Krt5")
An interactive version of the volcano plot above that includes the raw per sample values in a separate panel is possible via the glXYPlot function in the Glimma package.
library(Glimma)
group2 <- group
levels(group2) <- c("basal.lactate","basal.preg","basal.virgin","lum.lactate", "lum.preg", "lum.virgin")
glXYPlot(x=results$logFC, y=-log10(results$FDR),
xlab="logFC", ylab="B", main="B.PregVsLac",
counts=normCounts, groups=group2, status=de,
anno=ann, id.column="ENTREZID", folder="volcano")
This function creates an html page (./volcano/XY-Plot.html) with a volcano plot on the left and a plot showing the log-CPM per sample for a selected gene on the right. A search bar is available to search for genes of interest.
It might seem natural to add genomic locations to our annotation table, and possibly a bit odd that the org.Mm.eg.db package does not supply such mappings. In fact, there is a whole suite of package for performing this, and more-advanced queries. These are listed on the Bioconductor annotation page and have the prefix TxDb.
The package we will be using is TxDb.Mmusculus.UCSC.mm10.knownGene. Packages are available for other organisms and genome builds. It is even possible to build your own database if one does not exist. See vignette("GenomicFeatures") for details
source("http://www.bioconductor.org/biocLite.R")
biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")
## For Humans
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
We load the library in the usual fashion and create a new object to save some typing. As with the org. packages, we can query what columns are available with columns,
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
Loading required package: GenomicFeatures
Loading required package: GenomeInfoDb
Loading required package: GenomicRanges
tx <- TxDb.Mmusculus.UCSC.mm10.knownGene
columns(tx)
[1] "CDSCHROM" "CDSEND" "CDSID" "CDSNAME" "CDSSTART" "CDSSTRAND"
[7] "EXONCHROM" "EXONEND" "EXONID" "EXONNAME" "EXONRANK" "EXONSTART"
[13] "EXONSTRAND" "GENEID" "TXCHROM" "TXEND" "TXID" "TXNAME"
[19] "TXSTART" "TXSTRAND" "TXTYPE"
The select function is used in the same manner as the org.Mm.eg.db packages.
Challenge
Use the TxDb.Mmusculus.UCSC.mm10.knownGene package to retrieve the exon coordinates for the genes
50916,110308,12293
One of the real strengths of the txdb.. packages is the ability of interface with GenomicRanges, which is the object type used throughout Bioconductor to manipulate Genomic Intervals.
These object types permit us to perform common operations on intervals such as overlapping and counting.
library(GenomicRanges)
my.range <-GRanges("chr1", IRanges(start=1000,end=2000))
my.range
GRanges object with 1 range and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 [1000, 2000] *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
Let’s try some dummy intervals with random chromsome assignments and start positions.
set.seed(05052017)
chrs <- sample(paste0("chr",1:19),100,replace=TRUE)
startPos <- sample(1:5e6,100,replace=TRUE)
my.rangesA <- GRanges(chrs, IRanges(start=startPos,width=1000))
my.rangesA
GRanges object with 100 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr17 [2732950, 2733949] *
[2] chr3 [2250749, 2251748] *
[3] chr15 [2581981, 2582980] *
[4] chr4 [3439502, 3440501] *
[5] chr15 [2851346, 2852345] *
... ... ... ...
[96] chr8 [ 35930, 36929] *
[97] chr13 [ 644034, 645033] *
[98] chr15 [1278939, 1279938] *
[99] chr2 [3348995, 3349994] *
[100] chr5 [4484841, 4485840] *
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
There are a number of useful functions for calculating properties of the data (such as coverage or sorting). It is even possible to convert between different chromosome naming conventions.
width(my.rangesA)
[1] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
[22] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
[43] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
[64] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
[85] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
sort(my.rangesA)
GRanges object with 100 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr17 [ 571363, 572362] *
[2] chr17 [1125242, 1126241] *
[3] chr17 [1338513, 1339512] *
[4] chr17 [2732950, 2733949] *
[5] chr3 [ 609141, 610140] *
... ... ... ...
[96] chr14 [1658130, 1659129] *
[97] chr14 [3037287, 3038286] *
[98] chr14 [3805150, 3806149] *
[99] chr14 [4853072, 4854071] *
[100] chr7 [3795787, 3796786] *
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
coverage(my.rangesA)
RleList of length 18
$chr17
integer-Rle of length 2733949 with 8 runs
Lengths: 571362 1000 552879 1000 212271 1000 1393437 1000
Values : 0 1 0 1 0 1 0 1
$chr3
integer-Rle of length 3387478 with 6 runs
Lengths: 609140 1000 1640608 1000 1134730 1000
Values : 0 1 0 1 0 1
$chr15
integer-Rle of length 4431853 with 20 runs
Lengths: 21661 1000 135948 1000 133119 1000 ... 210004 1000 184691 1000 390271 1000
Values : 0 1 0 1 0 1 ... 0 1 0 1 0 1
$chr4
integer-Rle of length 4744112 with 10 runs
Lengths: 271227 1000 1796879 1000 1101428 1000 266967 1000 1302611 1000
Values : 0 1 0 1 0 1 0 1 0 1
$chr6
integer-Rle of length 4526564 with 12 runs
Lengths: 1810833 1000 571823 1000 1073193 1000 554296 1000 141906 1000 368513 1000
Values : 0 1 0 1 0 1 0 1 0 1 0 1
...
<13 more elements>
seqlevelsStyle(my.rangesA)
[1] "UCSC"
keepSeqlevels(my.rangesA,"chr19")
GRanges object with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr19 [3107336, 3108335] *
[2] chr19 [3918859, 3919858] *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
#seqlevelsStyle(my.rangesA) <- "Ensembl"
#my.rangesA
my.rangesB <- shift(my.rangesA, shift = 500)
my.rangesB <- resize(my.rangesB,width =sample(1:1000,100,replace = TRUE))
my.rangesB
GRanges object with 100 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr17 [2733450, 2734037] *
[2] chr3 [2251249, 2251485] *
[3] chr15 [2582481, 2582952] *
[4] chr4 [3440002, 3440135] *
[5] chr15 [2851846, 2852015] *
... ... ... ...
[96] chr8 [ 36430, 37283] *
[97] chr13 [ 644534, 645469] *
[98] chr15 [1279439, 1280030] *
[99] chr2 [3349495, 3350125] *
[100] chr5 [4485341, 4486290] *
-------
seqinfo: 18 sequences from an unspecified genome; no seqlengths
findOverlaps(my.rangesA,my.rangesB)
Hits object with 100 hits and 0 metadata columns:
queryHits subjectHits
<integer> <integer>
[1] 1 1
[2] 2 2
[3] 3 3
[4] 4 4
[5] 5 5
... ... ...
[96] 96 96
[97] 97 97
[98] 98 98
[99] 99 99
[100] 100 100
-------
queryLength: 100 / subjectLength: 100
It is quite straightforward to translate the output of a select query into a GenomicFeatures object. However, several convenience functions exist to retrieve the structure of every gene for a given organism in one object.
The output of exonsBy is a list, where each item in the list is the exon co-ordinates of a particular gene.
exo <- exonsBy(tx,"gene")
exo
GRangesList object of length 24116:
$100009600
GRanges object with 7 ranges and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr9 [21062393, 21062717] - | 134539 <NA>
[2] chr9 [21062894, 21062987] - | 134540 <NA>
[3] chr9 [21063314, 21063396] - | 134541 <NA>
[4] chr9 [21066024, 21066377] - | 134542 <NA>
[5] chr9 [21066940, 21067925] - | 134543 <NA>
[6] chr9 [21068030, 21068117] - | 134544 <NA>
[7] chr9 [21073075, 21075496] - | 134546 <NA>
$100009609
GRanges object with 6 ranges and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
[1] chr7 [84940169, 84941088] - | 109989 <NA>
[2] chr7 [84943141, 84943264] - | 109990 <NA>
[3] chr7 [84943504, 84943722] - | 109991 <NA>
[4] chr7 [84946200, 84947000] - | 109992 <NA>
[5] chr7 [84947372, 84947651] - | 109993 <NA>
[6] chr7 [84963816, 84964009] - | 109994 <NA>
$100009614
GRanges object with 1 range and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
[1] chr10 [77711446, 77712009] + | 143986 <NA>
...
<24113 more elements>
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
To access the structure of a particular gene, we can use the [[ syntax with the name of the gene (Entrez gene ID) within quote marks.
exo[["12992"]]
GRanges object with 14 ranges and 2 metadata columns:
seqnames ranges strand | exon_id exon_name
<Rle> <IRanges> <Rle> | <integer> <character>
[1] chr5 [87808121, 87808165] + | 68233 <NA>
[2] chr5 [87809897, 87809959] + | 68234 <NA>
[3] chr5 [87813088, 87813114] + | 68235 <NA>
[4] chr5 [87813838, 87813861] + | 68236 <NA>
[5] chr5 [87813941, 87813982] + | 68237 <NA>
... ... ... ... . ... ...
[10] chr5 [87819066, 87819095] + | 68242 <NA>
[11] chr5 [87819066, 87819110] + | 68243 <NA>
[12] chr5 [87820951, 87820995] + | 68244 <NA>
[13] chr5 [87822212, 87822322] + | 68245 <NA>
[14] chr5 [87824139, 87824421] + | 68246 <NA>
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
It is also possible to save the results of a Bioconductor analysis in a browser to enable interactive analysis and integration with other data types, or sharing with collaborators. For instance, we might want a browser track to indicate where our differentially-expressed genes are located. We shall use the bed format to display these locations. We will annotate the ranges with information from our analysis such as the fold-change and significance.
First we create a data frame for just the DE genes.
sigGenes <- results.annotated[detags,]
sigGenes
At the moment, we have a GenomicFeatures object that represents every exon. However, we do not need this level of granularity for the bed output, so we will collapse to a single region for each gene. First we the range function to obtain a single range for every gene and tranform to a more convenient object with unlist.
range(exo)
GRangesList object of length 24116:
$100009600
GRanges object with 1 range and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr9 [21062393, 21075496] -
$100009609
GRanges object with 1 range and 0 metadata columns:
seqnames ranges strand
[1] chr7 [84940169, 84964009] -
$100009614
GRanges object with 1 range and 0 metadata columns:
seqnames ranges strand
[1] chr10 [77711446, 77712009] +
...
<24113 more elements>
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
exoRanges <- unlist(range(exo))
sigRegions <- exoRanges[na.omit(match(sigGenes$ENTREZID, names(exoRanges)))]
sigRegions
GRanges object with 4393 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
497097 chr1 [3214482, 3671498] -
20671 chr1 [4490928, 4497354] -
58175 chr1 [4909576, 5070285] -
76187 chr1 [9548046, 9577968] +
72481 chr1 [9560833, 9631092] -
... ... ... ...
195727 chrX [161836430, 162159441] -
108012 chrX [163909017, 163933666] +
56078 chrX [163976822, 164028010] -
54156 chrX [166523007, 166585716] -
333605 chrX [167471306, 168577233] -
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
Rather than just representing the genomic locations, the .bed format is also able to colour each range according to some property of the analysis (e.g. direction and magnitude of change) to help highlight particular regions of interest. A score can also be displayed when a particular region is clicked-on. A useful propery of GenomicRanges is that we can attach metadata to each range using the mcols function. The metadata can be supplied in the form of a data frame.
sigRegions
GRanges object with 4393 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
497097 chr1 [3214482, 3671498] -
20671 chr1 [4490928, 4497354] -
58175 chr1 [4909576, 5070285] -
76187 chr1 [9548046, 9577968] +
72481 chr1 [9560833, 9631092] -
... ... ... ...
195727 chrX [161836430, 162159441] -
108012 chrX [163909017, 163933666] +
56078 chrX [163976822, 164028010] -
54156 chrX [166523007, 166585716] -
333605 chrX [167471306, 168577233] -
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
mcols(sigRegions) <- sigGenes[match(names(sigRegions), rownames(sigGenes)),]
sigRegions
GRanges object with 4393 ranges and 8 metadata columns:
seqnames ranges strand | logFC logCPM LR
<Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric>
497097 chr1 [3214482, 3671498] - | -10.947240 2.5236515 23.590694
20671 chr1 [4490928, 4497354] - | -2.673131 1.2418640 10.587241
58175 chr1 [4909576, 5070285] - | 4.471434 1.1240115 14.373441
76187 chr1 [9548046, 9577968] + | 3.033003 2.4071013 8.475630
72481 chr1 [9560833, 9631092] - | 2.136618 -0.2472752 6.137581
... ... ... ... . ... ... ...
195727 chrX [161836430, 162159441] - | -4.692618 3.1077765 13.912120
108012 chrX [163909017, 163933666] + | -2.176996 2.3209830 10.201704
56078 chrX [163976822, 164028010] - | 2.588361 5.1331379 6.923093
54156 chrX [166523007, 166585716] - | -6.964350 4.3915759 19.444204
333605 chrX [167471306, 168577233] - | -3.896071 0.6561936 6.917891
PValue FDR ENTREZID SYMBOL
<numeric> <numeric> <character> <character>
497097 1.191624e-06 0.0004377961 497097 Xkr4
20671 1.138708e-03 0.0078553551 20671 Sox17
58175 1.499018e-04 0.0020196485 58175 Rgs20
76187 3.599356e-03 0.0182496716 76187 Adhfe1
72481 1.323382e-02 0.0468231045 72481 2610203C22Rik
... ... ... ... ...
195727 1.915593e-04 0.0023438662 195727 Nhs
108012 1.403109e-03 0.0091706934 108012 Ap1s2
56078 8.508968e-03 0.0338900535 56078 Car5b
54156 1.035817e-05 0.0005205134 54156 Egfl6
333605 8.533757e-03 0.0339270998 333605 Frmpd4
GENENAME
<character>
497097 X-linked Kx blood group related 4
20671 SRY (sex determining region Y)-box 17
58175 regulator of G-protein signaling 20
76187 alcohol dehydrogenase, iron containing, 1
72481 RIKEN cDNA 2610203C22 gene
... ...
195727 Nance-Horan syndrome (human)
108012 adaptor-related protein complex 1, sigma 2 subunit
56078 carbonic anhydrase 5b, mitochondrial
54156 EGF-like-domain, multiple 6
333605 FERM and PDZ domain containing 4
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
The metadata we have added can also by used as a means to interrogate the ranges; as if the data were contained in a data frame.
sigRegions[order(sigRegions$LR,decreasing = TRUE)]
GRanges object with 4393 ranges and 8 metadata columns:
seqnames ranges strand | logFC logCPM LR PValue
<Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric> <numeric>
110308 chr15 [101707070, 101712891] - | -8.940578 10.264297 24.89789 6.044844e-07
50916 chr13 [ 73260497, 73269620] + | -8.636503 5.749781 24.80037 6.358512e-07
12293 chr5 [ 15934691, 16374511] + | -8.362247 6.794788 24.68526 6.749827e-07
56069 chr18 [ 61687915, 61692537] + | -8.419433 6.124377 24.41532 7.764861e-07
24117 chr10 [121034004, 121100642] + | -9.290691 6.757163 24.32506 8.137331e-07
... ... ... ... . ... ... ... ...
75723 chr9 [ 14541967, 14643529] - | -1.373919 8.545179 5.988178 0.01440207
238384 chr12 [102129419, 102267091] + | -3.283021 -2.062221 5.986796 0.01441336
72543 chr2 [ 33729956, 33887946] - | -1.146311 5.243682 5.986777 0.01441351
268301 chr10 [ 59221922, 59226433] + | 1.393078 4.386359 5.984774 0.01442989
20931 chr2 [ 26916421, 26920170] + | -1.132158 4.526824 5.981059 0.01446032
FDR ENTREZID SYMBOL
<numeric> <character> <character>
110308 0.0004377961 110308 Krt5
50916 0.0004377961 50916 Irx4
12293 0.0004377961 12293 Cacna2d1
56069 0.0004377961 56069 Il17b
24117 0.0004377961 24117 Wif1
... ... ... ...
75723 0.04982713 75723 Amotl1
238384 0.04984490 238384 Slc24a4
72543 0.04984490 72543 Mvb12b
268301 0.04989062 268301 Sowahc
20931 0.04998488 20931 Surf2
GENENAME
<character>
110308 keratin 5
50916 Iroquois related homeobox 4 (Drosophila)
12293 calcium channel, voltage-dependent, alpha2/delta subunit 1
56069 interleukin 17B
24117 Wnt inhibitory factor 1
... ...
75723 angiomotin-like 1
238384 solute carrier family 24 (sodium/potassium/calcium exchanger), member 4
72543 multivesicular body subunit 12B
268301 sosondowah ankyrin repeat domain family member C
20931 surfeit gene 2
-------
seqinfo: 66 sequences (1 circular) from mm10 genome
For visualisation purposes, we are going to restrict the data to genes that are located on chromosomes 1 to 19 and the sex chromosomes. This can be done with the keepSeqLevels function.
seqlevels(sigRegions)
[1] "chr1" "chr2" "chr3"
[4] "chr4" "chr5" "chr6"
[7] "chr7" "chr8" "chr9"
[10] "chr10" "chr11" "chr12"
[13] "chr13" "chr14" "chr15"
[16] "chr16" "chr17" "chr18"
[19] "chr19" "chrX" "chrY"
[22] "chrM" "chr1_GL456210_random" "chr1_GL456211_random"
[25] "chr1_GL456212_random" "chr1_GL456213_random" "chr1_GL456221_random"
[28] "chr4_GL456216_random" "chr4_GL456350_random" "chr4_JH584292_random"
[31] "chr4_JH584293_random" "chr4_JH584294_random" "chr4_JH584295_random"
[34] "chr5_GL456354_random" "chr5_JH584296_random" "chr5_JH584297_random"
[37] "chr5_JH584298_random" "chr5_JH584299_random" "chr7_GL456219_random"
[40] "chrX_GL456233_random" "chrY_JH584300_random" "chrY_JH584301_random"
[43] "chrY_JH584302_random" "chrY_JH584303_random" "chrUn_GL456239"
[46] "chrUn_GL456359" "chrUn_GL456360" "chrUn_GL456366"
[49] "chrUn_GL456367" "chrUn_GL456368" "chrUn_GL456370"
[52] "chrUn_GL456372" "chrUn_GL456378" "chrUn_GL456379"
[55] "chrUn_GL456381" "chrUn_GL456382" "chrUn_GL456383"
[58] "chrUn_GL456385" "chrUn_GL456387" "chrUn_GL456389"
[61] "chrUn_GL456390" "chrUn_GL456392" "chrUn_GL456393"
[64] "chrUn_GL456394" "chrUn_GL456396" "chrUn_JH584304"
sigRegions <- keepSeqlevels(sigRegions, paste0("chr", c(1:19,"X","Y")))
We will now create a score from the p-values that will displayed under each region, and colour scheme for the regions based on the fold-change. For the score we can use the \(-log_{10}\) of the adjusted p-value as before
Score <- -log10(sigRegions$FDR)
colorRampPalette is a useful function in base R for constructing a palette between two extremes. When choosing colour palettes, make sure they are colour blind friendly. The red / green colour scheme traditionally-applied to microarrays is a bad choice.
We will also truncate the fold-changes to between -5 and 5 to and divide this range into 10 equal bins
rbPal <-colorRampPalette(c("red", "blue"))
logfc <- pmax(sigRegions$logFC, -5)
logfc <- pmin(logfc , 5)
Col <- rbPal(10)[as.numeric(cut(logfc, breaks = 10))]
The colours and score have to be saved in the GRanges object as score and itemRgb columns respectively, and will be used to construct the browser track. The rtracklayer package can be used to import and export browsers tracks.
Now we can export the signifcant results from the DE analysis as a .bed track using rtracklayer. You can load the resulting file in IGV, if you wish.
mcols(sigRegions)$score <- Score
mcols(sigRegions)$itemRgb <- Col
sigRegions
GRanges object with 4392 ranges and 10 metadata columns:
seqnames ranges strand | logFC logCPM LR
<Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric>
497097 chr1 [3214482, 3671498] - | -10.947240 2.5236515 23.590694
20671 chr1 [4490928, 4497354] - | -2.673131 1.2418640 10.587241
58175 chr1 [4909576, 5070285] - | 4.471434 1.1240115 14.373441
76187 chr1 [9548046, 9577968] + | 3.033003 2.4071013 8.475630
72481 chr1 [9560833, 9631092] - | 2.136618 -0.2472752 6.137581
... ... ... ... . ... ... ...
195727 chrX [161836430, 162159441] - | -4.692618 3.1077765 13.912120
108012 chrX [163909017, 163933666] + | -2.176996 2.3209830 10.201704
56078 chrX [163976822, 164028010] - | 2.588361 5.1331379 6.923093
54156 chrX [166523007, 166585716] - | -6.964350 4.3915759 19.444204
333605 chrX [167471306, 168577233] - | -3.896071 0.6561936 6.917891
PValue FDR ENTREZID SYMBOL
<numeric> <numeric> <character> <character>
497097 1.191624e-06 0.0004377961 497097 Xkr4
20671 1.138708e-03 0.0078553551 20671 Sox17
58175 1.499018e-04 0.0020196485 58175 Rgs20
76187 3.599356e-03 0.0182496716 76187 Adhfe1
72481 1.323382e-02 0.0468231045 72481 2610203C22Rik
... ... ... ... ...
195727 1.915593e-04 0.0023438662 195727 Nhs
108012 1.403109e-03 0.0091706934 108012 Ap1s2
56078 8.508968e-03 0.0338900535 56078 Car5b
54156 1.035817e-05 0.0005205134 54156 Egfl6
333605 8.533757e-03 0.0339270998 333605 Frmpd4
GENENAME score itemRgb
<character> <numeric> <character>
497097 X-linked Kx blood group related 4 3.358728 #FF0000
20671 SRY (sex determining region Y)-box 17 2.104834 #C60038
58175 regulator of G-protein signaling 20 2.694724 #0000FF
76187 alcohol dehydrogenase, iron containing, 1 1.738745 #1C00E2
72481 RIKEN cDNA 2610203C22 gene 1.329540 #3800C6
... ... ... ...
195727 Nance-Horan syndrome (human) 2.630067 #FF0000
108012 adaptor-related protein complex 1, sigma 2 subunit 2.037598 #C60038
56078 carbonic anhydrase 5b, mitochondrial 1.469928 #3800C6
54156 EGF-like-domain, multiple 6 3.283568 #FF0000
333605 FERM and PDZ domain containing 4 1.469453 #E2001C
-------
seqinfo: 21 sequences from mm10 genome
library(rtracklayer)
export(sigRegions , con = "topHits.bed")
As we have been using counts as our starting point, we haven’t investigated the aligned reads from our experiment, and how they are represented. As you may be aware, aligned reads are usually stored in a bam file that can be manipulated with open-source command-line tools such as samtools and picard. Bioconductor provide a low-level interface to bam/sam files in the form of the Rsamtools package. The GenomicAlignments package can also be used to retrieve the reads mapping to a particular genomic region in an efficient manner.
library(GenomicAlignments)
generegion <- exo[["50916"]]
getwd()
bam <- readGAlignments(file="bam/MCL1.DG.bam",
param=ScanBamParam(which=generegion))
bam
There are various visualisation options for aligned reads and range data. We will use the ggbio package, which first requires some discussion of the ggplot2 plotting package.
The ggplot2 package has emerged as an attractive alternative to the traditional plots provided by base R. A full overview of all capabilities of the package is available from the cheatsheet.
A simple scatter plot, equivalent to plotSmear from before, can be generated as follows:-
library(ggplot2)
ggplot(results, aes(x = logCPM, y=logFC)) + geom_point()
In brief:-
results is our data frame containing the variables we wish to plotaes creates a mpping between the variables in our data frame to the aesthetic proprties of the plot
logCPM, y-axis is mapped to logFCgeom_point specifies the particular type of plot we want (in this case a scatter plot)
The real advantage of ggplot2 is the ability to change the appearance of our plot by mapping other variables to aspects of the plot. For example, we could colour the points based on a p-value cut-off. The colours are automatically chosen by ggplot2, but we can specifiy particular values.
ggplot(results, aes(x = logCPM, y=logFC,col=FDR < 0.05)) + geom_point()
ggplot(results, aes(x = logCPM, y=logFC,col=FDR < 0.05)) + geom_point(alpha=0.4) + scale_colour_manual(values=c("black","red"))
ggplot(results, aes(x = logCPM, y=logFC,col=FDR < 0.05)) + geom_point() + scale_colour_manual(values=c("black","red"))
The volcano plot can be constructed in a similar manner
ggplot(results, aes(x = logFC, y=-log10(FDR))) + geom_point()
We will now take a brief look at one of the visualisation packages in Bioconductor that takes advantage of the GenomicRanges and GenomicFeatures object-types. In this section we will show a worked example of how to combine several types of genomic data on the same plot. The documentation for ggbio is very extensive and contains lots of examples.
http://www.tengfei.name/ggbio/docs/
The Gviz package is another Bioconductor package that specialising in genomic visualisations, but we will not explore this package in the course.
The Manhattan plot is a common way of visualising genome-wide results, and this is implemented as the plotGrandLinear function. We have to supply a value to display on the y-axis using the aes function, which is inherited from ggplot2. The positioning of points on the x-axis is handled automatically by ggbio, using the ranges information to get the genomic coordinates of the ranges of interest.
To stop the plots from being too cluttered we will consider the top 200 genes only.
library(ggbio)
Need specific help about ggbio? try mailing
the maintainer or visit http://tengfei.github.com/ggbio/
Attaching package: 'ggbio'
The following objects are masked from 'package:ggplot2':
geom_bar, geom_rect, geom_segment, ggsave, stat_bin, stat_identity, xlim
top200 <- sigRegions[order(sigRegions$LR,decreasing = TRUE)[1:200]]
plotGrandLinear(top200 , aes(y = logFC))
using coord:genome to parse x scale
ggbio has alternated the colours of the chromosomes. However, an appealing feature of ggplot2 is the ability to map properties of your plot to variables present in your data. For example, we could create a variable to distinguish between up- and down-regulated genes. The variables used for aesthetic mapping must be present in the mcols section of your ranges object.
mcols(top200)$UpRegulated <- mcols(top200)$logFC > 0
plotGrandLinear(top200, aes(y = logFC, col = UpRegulated))
using coord:genome to parse x scale
plotGrandLinear is a special function in ggbio with preset options for the manhattan style of plot. More often, users will call the autoplot function and ggbio will choose the most appropriate layout. One such layout is the karyogram.
autoplot(top200,layout="karyogram",aes(color=UpRegulated,
fill=UpRegulated))
Scale for 'x' is already present. Adding another scale for 'x', which will replace the
existing scale.
Scale for 'x' is already present. Adding another scale for 'x', which will replace the
existing scale.
ggbio is also able to plot the structure of genes according to a particular model represented by a GenomicFeatures object.
autoplot(tx, which=exo[["24117"]])
Parsing transcripts...
Parsing exons...
Parsing cds...
Parsing utrs...
------exons...
------cdss...
------introns...
------utr...
aggregating...
Done
Constructing graphics...
We can even plot the location of sequencing reads if they have been imported using readGAlignments function (or similar).
myreg <- flank(reduce(exo[["24117"]]), 1000, both = T)
bam <- readGAlignments(file="bam/MCL1.DG.bam",
param=ScanBamParam(which=myreg),use.names = TRUE)
autoplot(bam,which=myreg)
autoplot(bam , stat = "coverage")
Like ggplot2, ggbio plots can be saved as objects that can later be modified, or combined together to form more complicated plots. If saved in this way, the plot will only be displayed on a plotting device when we query the object. This strategy is useful when we want to add a common element (such as an ideogram) to a plot composition and don’t want to repeat the code to generate the plot every time.
#idPlot <- plotIdeogram(genome = "mm10",subchr = "chr1")
#idPlot
geneMod <- autoplot(tx, which = myreg)
reads1 <- autoplot(bam, stat = "coverage")
tracks(geneMod ,reads1)